using System;
using System.Collections.Generic;
using System.Linq;

using Atomic.Libraries;
using Atomic.Libraries.Mathematics;
using Atomic.Libraries.Mathematics.Integration;

namespace Atomic.Thermodynamics.Electrons
{
	/// <summary>
	/// Computes the electronic free energy. The non-excited ground state energy subtracted from the free energy.
	/// </summary>
	public class ElectronicPotential : IPotential
	{
		private IElectronicStateDensityPoint[] densities;
		private double fermiEnergy, electrons;
		[NonSerialized] private InnerEvaluator evaluator;
		[NonSerialized] private RollingCache<double, Tuple<double, double, double>> values;

		public ElectronicPotential(IElectronicStateDensity stateDensity)
			: this(stateDensity, false)
		{
		}

		public ElectronicPotential(IElectronicStateDensity stateDensity, bool calculateFermiEnergy)
		{
			//double[] energies, densities0, densities1; // energies, and linear interpolation coefficients
			densities = stateDensity.Points.ToArray();
			fermiEnergy = stateDensity.FermiEnergy;
			electrons = stateDensity.Electrons;

			if (calculateFermiEnergy)
			{
				// Recomputes the Fermi energy such that the number of electrons is consistent.
				fermiEnergy = new InnerEvaluator(densities, fermiEnergy, electrons).CalculateFermiEnergy();
			}
		}

		public double FreeEnergy(double temperature)
		{
			// Compute for each temperature. No linear interpolation! Plotting is very easy now, e.g. t => p.FreeEnergy(t).
			return Compute(temperature).Item1;
		}

		public double Entropy(double temperature)
		{
			return Compute(temperature).Item2;
		}

		public double HeatCapacity(double temperature)
		{
			return Compute(temperature).Item3;
		}

		/// <summary>
		/// The interpolated density of states function used for the calculations.
		/// </summary>
		public double StateDensity(double energy)
		{
			if (energy < densities[0].Energy)
			{
				return 0.0;
			}

			for (int i = 1; i < densities.Length; i++)
			{
				double e0 = densities[i - 1].Energy;
				double e = densities[i].Energy;

				// Constant term and slope of linear interpolation.
				double d = (densities[i].IntegratedDensity - densities[i - 1].IntegratedDensity) / (e - e0);

				if (energy < e)
				{
					return d;
				}
			}

			// Beyond last point.
			return 0.0;
		}

		public double Electrons()
		{
			return new InnerEvaluator(densities, fermiEnergy, electrons).AverageElectronCount(0.0, double.PositiveInfinity);
		}

		private Tuple<double, double, double> Compute(double temperature)
		{
			if (values == null)
			{
				values = new RollingCache<double, Tuple<double, double, double>>(10000);
			}

			if (evaluator == null)
			{
				evaluator = new InnerEvaluator(densities, fermiEnergy, electrons);
			}

			Tuple<double, double, double> value;
			if (!values.TryGetValue(temperature, out value))
			{
				double f, s, cv;
				evaluator.Compute(temperature, out f, out s, out cv);

				// Cache computed free energy, entropy, and heat capacity as a tuple.
				value = Tuple.Create<double, double, double>(f, s, cv);
				values.Add(temperature, value);
			}

			return value;
		}

		[Serializable]
		private class InnerEvaluator
		{
			private static GaussLegendreQuadrature quadrature;
			private IElectronicStateDensityPoint[] densities;
			private double fermiEnergy, electrons;

			static InnerEvaluator()
			{
				// Currently just using numerical integration. Prepare quadrature once.
				quadrature = new GaussLegendreQuadrature(5);
			}

			public InnerEvaluator(IElectronicStateDensityPoint[] densities, double fermiEnergy, double electrons)
			{
				this.densities = densities;
				this.fermiEnergy = fermiEnergy;
				this.electrons = electrons;
			}

			public void Compute(double t, out double f, out double s, out double cv)
			{
				if (t == 0.0)
				{
					// Small deviations in the exact and approximate formulas at zero-temperature may exist. Return exact values.
					f = 0.0;
					s = double.NaN; // ??
					cv = double.NaN; // ??
					return;
				}

				double kb = 8.617332478e-5; // eV/K
				double beta = 1.0 / (kb * t);

				// Determine chemical potential required to maintain constant number of electrons (we don't require constant spin-polarization).
				double mu = ElectronicChemicalPotential(beta);

				// Subdivisions.
				int m = 15;
				double delta = 0.5 / beta;

				// Internal energy, internal energy at zero temperature, and entropy.
				double u = 0.0;
				double u0 = 0.0;
				s = 0.0;

				double e0 = densities[0].Energy - fermiEnergy;
				for (int i = 1; i < densities.Length; i++)
				{
					double e = densities[i].Energy - fermiEnergy;

					// The integrated density of state function is linearly interpolated. Hence the density itself is piece-wise constant.
					double d = (densities[i].IntegratedDensity - densities[i - 1].IntegratedDensity) / (e - e0);

					if (d == 0.0)
					{
						// Don't waste time on this. Happens often.
					}
					else
					{
						// Introduce a few more steps around the chemical potential (if overlapping with current interval).
						// Makes the code more messy, but necessary to obtain numerically stable estimates.
						double ej0 = e0;
						for (int j = -m; j <= m + 1; j++)
						{
							// Never go below original interval.
							double ej = Math.Max(e0, mu + j * delta);

							if (j > m || ej > e)
							{
								// Always end at original interval. But don't go beyond.
								ej = e;
							}

							if (ej != ej0)
							{
								// Here comes the actual integration.

								// Internal energy $U=\int d\eps n_F(\eps)\eps$ and grand-canonical free energy $F=-\beta^{-1}...$.
								u += quadrature.Integrate(x => d * x * FermiDiracDistribution(x, mu, beta), ej0, ej);
								s += quadrature.Integrate(x => { double ff = FermiDiracDistribution(x, mu, beta); return d * ((ff != 0.0 ? (ff * Math.Log(ff)) : 0.0) + (ff != 1.0 ? ((1.0 - ff) * Math.Log(1.0 - ff)) : 0.0)); }, ej0, ej);
								//cv += density.Value(e) * Math.Exp((e - mu) / (kb * t)) * Math.Pow(ff * (e - mu) / t, 2.0) / kb * dE; // from T*dS/dT
							}

							ej0 = ej;
						}
					}

					// Finally, compute internal energy at zero-temperature. Just the step function. Solve trivial integral in closed form.
					u0 += 0.5 * d * (Math.Min(e, 0.0) * Math.Min(e, 0.0) - Math.Min(e0, 0.0) * Math.Min(e0, 0.0));

					e0 = e;
				}

				u -= u0;
				s *= -kb;
				f = u - t * s;
				cv = double.NaN;
			}

			public double CalculateFermiEnergy()
			{
				// Recomputes the Fermi energy such that the number of electrons is consistent.
				return fermiEnergy + RootFinder.Bisection(mu => AverageElectronCount(mu, double.PositiveInfinity) - electrons, -0.6, 0.9);
			}

			public double ElectronicChemicalPotential(double beta)
			{
				// Determine number of electrons by integrating the density of states up to the Fermi energy
				// at zero temperature.
				double n = AverageElectronCount(0.0, double.PositiveInfinity);

				// Determine the chemical potential required to maintain this number of electrons. The search
				// interval is assumed to be wide enough, as the chemical potential doesn't vary much (otherwise
				// NaN is returned by the root finder algorithm).
				return RootFinder.Bisection(mu => AverageElectronCount(mu, beta) - n, -0.6, 0.9);
			}

			public double AverageElectronCount(double mu, double beta)
			{
				// Subdivisions.
				int m = 6;
				double delta = 1.0 / beta;
	
				// Number of particles.
				double n = 0.0;

				double e0 = densities[0].Energy - fermiEnergy;
				for (int i = 1; i < densities.Length; i++)
				{
					double e = densities[i].Energy - fermiEnergy;

					// The integrated density of state function is linearly interpolated. Hence the density itself is piece-wise constant.
					double d = (densities[i].IntegratedDensity - densities[i - 1].IntegratedDensity) / (e - e0);

					if (d == 0.0)
					{
						// Don't waste time on this. Happens often.
					}
					else if (double.IsPositiveInfinity(beta))
					{
						// Just the step function. Solve trivial integral in closed form.
						n += d * (Math.Min(e, mu) - Math.Min(e0, mu));
					}
					else
					{
						// Introduce a few more steps around the chemical potential (if overlapping with current interval).
						// Makes the code more messy, but necessary to obtain numerically stable estimates.
						double ej0 = e0;
						for (int j = -m; j <= m + 1; j++)
						{
							// Never go below original interval.
							double ej = Math.Max(e0, mu + j * delta);

							if (j > m || ej > e)
							{
								// Always end at original interval. But don't go beyond.
								ej = e;
							}

							if (ej != ej0)
							{
								// Here comes the actual integration.
								n += quadrature.Integrate(x => d * FermiDiracDistribution(x, mu, beta), ej0, ej);
							}

							ej0 = ej;
						}
					}

					e0 = e;
				}
			
				return n;
			}

			private static double FermiDiracDistribution(double e, double mu, double beta)
			{
				return 1.0 / (Math.Exp(beta * (e - mu)) + 1.0);
			}

			/*/// <summary>
			/// Incomplete Fermi-Dirac integral $F_0(x,b)=\int_b^\infty 1/(exp(t-x)+1)dt$.
			/// </summary>
			private static double IncompleteFermiDiracIntegral0(double x, double b)
			{
				//return Gsl.gsl_sf_fermi_dirac_inc_0(x, b);
				//return Gsl.gsl_sf_fermi_dirac_0(x - b);
				return Math.Log(1.0 + Math.Exp(x - b));
			}

			/// <summary>
			/// Incomplete Fermi-Dirac integral $F_0(x,a)-F_0(x,b)=\int_a^b 1/(exp(t-x)+1)dt$.
			/// </summary>
			private static double IncompleteFermiDiracIntegral0(double x, double a, double b)
			{
				// Determine common factor to divide out. Avoid this if both are negative.
				double x0 = Math.Max(0.0, Math.Max(x - a, x - b));
				return Math.Log((Math.Exp(-x0) + Math.Exp(x - a - x0)) / (Math.Exp(-x0) + Math.Exp(x - b - x0)));

				//return IncompleteFermiDiracIntegral0(x, a) - IncompleteFermiDiracIntegral0(x, b);
				//return Math.Log(1.0 + Math.Exp(x - a)) - Math.Log(1.0 + Math.Exp(x - b));
			}

			private static double IncompleteFermiDiracIntegral0(double x, double a, double b, double beta)
			{
				return IncompleteFermiDiracIntegral0(beta * x, beta * a, beta * b) / beta;
			}

			private static double IncompleteFermiDiracIntegral0_numerical(double x, double a, double b)
			{
				return GaussLegendreQuadrature.Integrate(t => 1.0 / (Math.Exp(t - x) + 1.0), a, b, 100);
			}

			/// <summary>
			/// Incomplete Fermi-Dirac integral $F_1(x,b)=\int_b^\infty t/(exp(t-x)+1)dt$.
			/// </summary>
			private static double IncompleteFermiDiracIntegral1(double x, double b)
			{
				// http://math.stackexchange.com/questions/471066/incomplete-fermi-dirac-integrals-and-polylogs
				double y = Math.Exp(b - x);
				//Console.WriteLine(Math.PI / 6.0 - b * b / 2.0 + x * x / 2.0 + b * Math.Log(1.0 + y) + SpecialFunctions.PolyLog(2, -y));
				//Console.WriteLine(Gsl.gsl_sf_fermi_dirac_1(x - b) + b * Gsl.gsl_sf_fermi_dirac_0(x - b));
				return Math.PI / 6.0 - b * b / 2.0 + x * x / 2.0 + b * Math.Log(1.0 + y) + SpecialFunctions.PolyLog(2, -y);
			}

			/// <summary>
			/// Incomplete Fermi-Dirac integral $F_1(x,a)-F_1(x,b)=\int_a^b t/(exp(t-x)+1)dt$.
			/// </summary>
			private static double IncompleteFermiDiracIntegral1(double x, double a, double b)
			{
				//return IncompleteFermiDiracIntegral1(x, a) - IncompleteFermiDiracIntegral1(x, b);
				return IncompleteFermiDiracIntegral1_numerical(x, a, b);
			}

			private static double IncompleteFermiDiracIntegral1(double x, double a, double b, double beta)
			{
				return IncompleteFermiDiracIntegral1(beta * x, beta * a, beta * b) / (beta * beta);
			}

			private static double IncompleteFermiDiracIntegral1_numerical(double x, double a, double b)
			{
				return GaussLegendreQuadrature.Integrate(t => t / (Math.Exp(t - x) + 1.0), a, b, 100);
			}

			private static double IncompleteFermiDiracIntegral2(double x, double a, double b, double beta)
			{
				return IncompleteFermiDiracIntegral2_numerical(beta * x, beta * a, beta * b) / (beta * beta * beta);
			}

			private static double IncompleteFermiDiracIntegral2_numerical(double x, double a, double b)
			{
				return GaussLegendreQuadrature.Integrate(t => t * t / (Math.Exp(t - x) + 1.0), a, b, 100);
			}*/
		}
	}
}
